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Abstract 



In this paper the unconditional stability of four well-known ADI schemes is analyzed in the 
application to time-dependent multidimensional diffusion equations with mixed derivative 
' terms. Necessary and sufficient conditions on the parameter 8 of each scheme are obtained 

that take into account the actual size of the mixed derivative coefficients. Our results gener- 
alize results obtained previously by Craig & Sneyd (1988) and In 't Hout & Welfert (2009). 
Numerical experiments are presented illustrating our main theorems. 

1 Introduction 

>- ! 

^ | 1.1 Multidimensional diffusion equations with mixed derivative terms 



This paper is concerned with stability in the numerical solution of initial-boundary value problems 
for time-dependent multidimensional diffusion equations containing mixed spatial-derivative terms, 



. du ^ 
O ■ "37 = dijUxiXj + d n u x lXl + a 2 2U X2X2 H h d kk u XkXk (1.1) 

(N ■ 61 i& 

Here k > 2 is any given integer and D — (<Ay)i<ij<fc denotes any given real, symmetric, positive 
semidefinite k x k matrix. The spatial domain is taken as ft — (0, l) fc and t > 0. 
. Multidimensional diffusion equations with mixed derivative terms play an important role, no- 

tably, in the field of financial option pricing theory. Here these terms represent correlations 
between underlying stochastic processes, such as for asset prices, volatilities or interest rates, see 
e.g. P[16l[T3[T8]. 

In this paper we consider for given 76 [0,1] the following useful condition on the relative size 
of the mixed derivative coefficients, 



\dij \ < 7 \J dudjj whenever 1 < i, j < k, i 7^ j, (1-2) 

Clearly, 7 = yields the smallest possible mixed derivative coefficients; they all vanish in this case. 
Subsequently, 7 = 1 admits the largest possible mixed derivative coefficients, given that matrix 
D is positive semidefinite. In most applications in financial mathematics it turns out that (|1.2I) 
holds with certain known < 7 < 1. More precisely, such 7 can be determined from the pertinent 
correlation coefficients. The aim of our present paper is to effectively use this knowledge about 7 
in the stability analysis of numerical schemes, discussed below, so as to arrive at weaker sufficient 
stability conditions than those obtained when considering 7 = 1. 
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1.2 Alternating Direction Implicit schemes 

Semidiscretization of initial-boundary value problems for multidimensional diffusion equations 
(jTTTJ) leads to initial value problems for very large systems of stiff ordinary differential equations, 

U'{t)=AU(t)+g(t) (t>0), U(0) = U . (1.3) 

Here A is a given real m x m matrix and g(t) (for t > 0) and Uq are given real raxl vectors. For 
the efficient numerical solution of (|1.3p . splitting schemes of the Alternating Direction Implicit 
(ADI) type form an attractive means. These schemes employ a splitting of the matrix A that 
corresponds to the different spatial directions, leading to a substantial reduction in the amount of 
computational work per time step compared to common implicit methods such as Crank-Nicolson. 
ADI schemes constitute a popular class of numerical methods for multidimensional equations in 
financial mathematics, where mixed derivative terms are ubiquitous. 

Originally ADI schemes were proposed and analyzed for equations without mixed derivatives, 
see e.g. [UI31[S]. To formulate ADI schemes in the presence of mixed derivative terms, consider a 
splitting of the matrix A into 

A = A +A 1 + --- + A k , (1.4) 

where Aq represents all mixed derivative terms in (|1.1[) and Aj represents the derivative term in 
the j-th spatial direction for j = 1, 2, . . . , fc. Split the vector function g in an analogous way, 

9 = 9o + 9i H \-9k, 

and define 

F(t, £) = AC + g(t) and Fj(t, £) = Aj£ + gj(t) whenever t > 0, f <E R m , < j < k. 

Let parameter 9 > be given. Let time step At > and temporal grid points t n = nAt with 
integer n > 0. For the numerical solution of the semidiscrete system (|1.3j) . we study in this paper 
four ADI schemes, each successively generating for n— 1, 2,3, . . . approximations U n to U(t n ): 

Douglas (Do) scheme 

' r = £/ n _i+AtF(*»_i,Z7„_i), 

< Y^Yj-i+OAtiFjit^Y^-Fjitn^Un-i)), j = 1, 2, . . . , k, (1.5) 

Craig-Sneyd ( CS) scheme 

' Y Q = U n -x + AtF{t n ^ 1 ,U n - 1 ), 
Yj = Yj^ + 9 At {Fj{t n , Yj) - f^tn-i, I7„_i)) , j = 1, 2, . . . , fc, 

< F = r + |At (F (t n ,y fc )-Fo(t„-i, Un-i)), (1.6) 
$5 = + #Ai (Fj (t„ , Yj ) - Fj-ftn-i, [/„_!)), j = 1, 2, . . . , fc, 

. c/„ = r fe 
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Modified Craig-Sneyd (MCS) scheme 

' K = £/n-i+A£F(t B _i, 





= Yj-i 


+ 9 At (F 3 {t n ,Y 3 ) - Fjitn-!, f/ n _i)) , j = 1, 2, . 


. . , fc, 


% 


= Y Q + 


9 At (F (t n , Y k ) - F (t„_i, , 




% 


= % + 


(i - 9) At {F(t n , Y k ) - F(i„_ l5 l/ n _i)) , 






= Yj-i 


+ 0Ai (F 3 -(t n , Fj) - Fj(t n ^i, U n -i)), j = 1,2,. 




u„ 


= % 







Hundsdorjer-Verwer (HV) scheme 

' Y Q = [/„_! + AtF(t n _i, C/„_i), 

^ = + {Fj(t n , Yj) - fXtn-i, C/„-i)) , j = 1, 2, . . . , fc, 

< ? =Fo + iAt(F(t„,n)-F(t„_ lj [/„_!)), (1.8) 

^ = + 0Ai (Fj - Fj (t ns y fc )), j = 1, 2, . . . , k, 

. U n = Y k . 

A study of the above four ADI schemes shows that the Aq part, representing all mixed derivative 
terms, is always treated in an explicit fashion. In line with the original ADI idea, the Aj parts for 
j = 1, 2, . . . , k are successively treated in an implicit fashion. 

In its general form (11.51) the Douglas scheme has been considered e.g. in [TU III] . Special 
instances include the well-known ADI schemes of Douglas & Rachford [5], where Fq — and 
9 = 1, and of Brian [5] and Douglas g], where F = and 9 = ±. McKee & Mitchell [TS] first 
proposed the Do scheme for equations (|1.1[) where a mixed derivative term is present .with k = 2. 

It is readily shown that if A is nonzero then the classical order of consistency of the Do 
scheme is just one, for any given 9. 

The CS, MCS and HV schemes can be viewed as different extensions to the Do scheme. They 
each perform a second explicit prediction stage followed by k implicit correction stages. The CS 
scheme was introduced in [3] and attains a classical order equal to two (independently of Ao) 
if 9 = T). The MCS scheme was defined in [TT] and possesses a classical order equal to two for 
any given 9. Note that the CS and MCS schemes are identical if 9 = h. The HV scheme was 
constructed in [13] and was proposed for equations with mixed derivative terms in [10 . Like the 
MCS scheme, the HV scheme possesses the favorable property of having a classical order equal to 
two for any given 9. 

For the stability analysis in this paper the following linear scalar test equation is relevant, 

U'{t) = (\ + \ 1 + --- + \ k )U(t), (1.9) 

with complex constants Aj (0 < j < k). Application of any given ADI scheme in the case of test 
equation (|1.9|) gives rise to a linear iteration of the form 

U„ = M(z , zx, . . . , z k ) C/ n _i . 

1 That is, the order for fixed, nonstiff ODEs. 
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Here M is a given, fixed, multivariate rational function and Zj = At Xj (0 < j < k). The iteration 
is stable if 

\M(z ,z u ...,z k )\<l. (1.10) 

Write 

Z= z 1 +z 2 + --- + z k and p=(l-6zi)(l-9z 2 )---(l-9z k ). (1.11) 



For the schemes (|L5]> . (|L6|) . (fTT7| . (|L8|I it is readily verified that M = R, S, S, T, respectively, 
where 

R(z , Zl ,...,z k ) = l + ^±f ) (1.12) 

P 



5(zo,zi, . . . ,z k ) = H h- j ' i 1 ' 13 ) 



^ + Z 1 £o(fO_+ z) 

p 2 p 2 

a , ^ i , z o + z , fl z o(z + z) n , (z + z) 2 
5(zo,^i, •••,%) = H 1-0 5 ^\2~°) 2 — 



9 1 V 2 / 9 



(1.14) 



s 2n + 2 Zn + Z 1 (zn + z) 2 

T(z ,z 1 ,...,z k ) = 1 + 2^ —2 ■ !- 15 

1.3 Finite difference discretization 

For the semidiscretization of (jl.ip we replace all spatial derivatives by second-order central finite 
differences on a rectangular grid with constant mesh width Ax t > in the Xi-direction (1 < i < k): 



(«...,) 



U£ +£i — 2ii£ + ui- 

(1 + Pij){ui +ei+ej + Uf-e.-eJ ~ (1 ~ Pij ) (W- ei +ej + U£+e,- ej ) 



{u XiXi ) e « - x -y 2 , (1.16a) 



1 ' ^AxiAxj 

/3ij(4:U£ - 2(ut +ei + Ui +ej +U£- ei + t^_ e ,)) . . . , 
+ 4A^ ' (L16b) 



Here t = (tx,£2, ■ ■ ■ ,£ k ) and the unit vectors e±,e2, ■ ■ ■ ,e k denote multi-indices, ui stands for 
u{l% Ax\, £2 Ax 2 , . . . ,£ k Ax k ,t) and ftj (1 < i ^ j < k) are given real parameters with = (iji. 
The righthand side of (|1.16b ) is the most general second-order finite difference formula for the 
mixed derivative u XiX . that is based on a centered 9-point stencil. If /3y = 0, then (|1.16b ) reduces 
to the standard 4-point formula 

/ \ Ue+ei+ej + Uf- e e — U(- e +e . — U( +e e 



****** )t~ AAxiAxj 
In the literature also the choices (3ij = — 1 and /3y = 1 are frequently considered. 



1.4 Stability analysis and outline of this paper 

The aim of this paper is to investigate the stability of the four ADI schemes (|1.5[) - (|1.8p in the 
application to semidiscretized multidimensional diffusion equations (jl.ll) where the condition p. 21) 
on the size of the mixed derivative coefficients is effectively taken into account. Our stability 
analysis is equivalent to the well-known von Neumann (Fourier) analysis. Accordingly, (11.11) 
is considered with constant diffusion matrix D and periodic boundary condition and stability 
pertains to the ^-norm. The obtained semidiscrete matrices Aj (0 < j < k) are then all normal 
and commuting and can therefore be unitarily diagonalized by a single matrix. For any given ADI 
scheme, one thus arrives at the stability requirement (|1 . 10|) with Aj eigenvalues of the matrices 
Aj(0<j<k). 
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Upon substituting discrete Fourier modes into (I1.16[) . it is readily shown that the scaled eigen- 
values Zj are given by 



z = yVjj djj [- sintfo sin^j + - cos^)(l - cos^-)] , (1.17a) 

Zj = — 2rjjdjj(l — cos</>j), j = 1, 2, . . . , k. (1.17b) 

The angles <pj are integer multiples of 2tt jrrij with rrij the dimension of the grid in the a; j -direction 
(1 < j < k). Further, 

At 

Tii = — A : — for 1 < i, j < k. 

J AxiAxj -in- 
stability results for ADI schemes pertinent to (|1.2[) in the case of 7 = 1 were derived by 
Craig & Sneyd [3] for the Do and CS schemes and by In 't Hout & Welfert pi] for the MCS and 
HV schemes. For any given scheme, both necessary and sufficient conditions on its parameter 9 
were obtained such that the stability requirement (| 1 . 10[) with (zq, Z\, . . . , Zk) given by (|1.17[) is 
unconditionally fulfilled, i.e., for all At > and all Axi > (1 < i < k). 

The stability analysis of ADI schemes based on (|1.2|l with arbitrary given 7 € [0, 1] was recently 
started in In 't Hout & Mishra [8]. For the MCS scheme and k = 2, the useful result was proved 
here that the lower bound 9 > max{|, 5(7 + 1)} is sufficient for unconditional stability (whenever 
0<7<1). 

The present paper substantially extends the work of 3, 8||TT] reviewed above. Section [2] con- 
tains the two main results of this paper. The first main result is Theorem 12 .31 which provides 
for each of the Do, CS, MCS, HV schemes in k = 2 and k = 3 spatial dimensions a sufficient 
condition on the parameter 9 for unconditional stability under (|1.2|) for arbitrary given 7 G [0, 1]. 
The second main result is Theorem 12.41 This yields for any of the four ADI schemes, any spatial 
dimension k > 2 and any 7 G [0, 1] a necessary condition on 9 for unconditional stability under 
(|1.2[) . For each scheme, the obtained necessary and sufficient conditions coincide if k — 2 or k = 3. 
Section [3] presents numerical illustrations to the two main theorems. The final Section 0] gives 
conclusions and issues for future research. 

2 Main results 

In the following we always make the minor assumption that the matrix B = (— /%)i<ij<fc with 
Pa = — 1 (1 < i < k) is positive semidefinite. Thus, in particular, \(3ij\ < 1 for all i,j. We note 
that this assumption on B is weaker than the corresponding assumption that was made in [llj . 

2.1 Preliminaries 

This section gives two lemmas that shall be used in the proofs of the main results below. 
Lemma 2.1 Let a, 8 be given real numbers with < 5 < 4. Consider the polynomial 
P(u 7 v 7 w) — a + u 2 + v + w 2 + uvw — 5(u + v + w) for u,v,w € R. 



Then 

if and only if 



P{u 1 v, w) > whenever u>0,v>0,w>0 



(<5 + l)(3-2V^+l) > 1-a and 5 2 < 2a. (2.1) 
Proof The critical points (u, v, w) of P are given by the equations 

2u + vw = 2v + uw = 2w + uv = 5. 
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A straightforward analysis, using < S < 4, shows that there is precisely one critical point in the 
domain u > 0, v > 0, w > 0, and it is given by u = v = w = s/S + 1 — 1 =: u* . Inserting this into 
P and rewriting yields 

P{u,u\u) = (u* + 1) 3 -3(6+l)u* +a-l 

= (VsTTy - 3 (8 + 1) (VsTT - 1^) +a-l 

= (5 + 1) (3 - 2V6 + 1) + a - 1. 

Hence P(u*,u*, u*) > if and only if the first inequality in (|2.ip holds. 

Consider next the polynomial P on the boundary of the pertinent domain. It is clear that 
P(u, v,w) > if u,v,w — > 00. Next, on the boundary part w > 0, v > 0, w — there holds 

P(u,v,0) = a + u 2 + v 2 - 6{u + v) 

= {u - \5f + {v - \8f + a - \5 2 . 

Thus P(u, v, 0) > (whenever u > 0, v > 0) if and only if the second inequality in (|2.1[) holds. By 
symmetry, the result for the other two boundary parts is the same, which completes the proof. 



The subsequent analysis relies upon four key properties of the scaled eigenvalues (|1.17JI . 
Lemma 2.2 Let zq, z\, . . . ,Zk be given by Let 7 € [0, 1] and assume hi. 2}) holds. Then: 



all Zj are real, (2.2a) 

Zj < for 1 < j < fc, (2.2b) 

z + z o <0, (2.2c) 

\z \<iJ2^^- ( 2 - 2d ) 



Proof Properties (|2.2h ). fl2.2b ) are obvious. We consider (|2.2b ) and (|2.2H ). Using r,j = ^/rurJJ 
and the simple identity 

(sin0) 2 + (1 - cos0) 2 = 2(1 - cos^) 

one readily verifies that 

fc 

— (z + zq) = rij dij [sin^j sin cf>j — (1 — cos (f>i)(l — cos (f>j)] = u T Du + v T (Do B)v, 
where o denotes the Hadamard product of two matrices and 

It is well-known that the Hadamard product of two positive semidefinite matrices is also positive 
semidefmite, see e.g. [7], and hence, z + zq < is obtained. 

Next, using the Cauchy-Schwarz inequality and \ < 1, it directly follows that 

I sin <j>i sin <j)j — ftj (1 — cos 4>i){l — cos (f>j)\ < y/2(l — cos <pi) ■ J 2(1 — cos (f>j). 
Together with (|1.2[) this yields 

I -2o I < 1^2 V^Vj V d ndj3 \/2(l - cos^) • ^2(1 - cos^) = 7^ ■ 



G 



In the following we shall use the notation yj — ^J~9zj (for 1 < j < k). Then 

k k 

p=l[(l+y])>l and z=--Y,y 2 j (2-3) 

3=1 3=1 

and by §FB& ) there holds 

z+z Q > z - |z i > Y^zj -7XJv^ = ~# X!-"- • " X!'/'-"' • ( 2 - 4 ) 

3=1 \3=1 / 

2.2 Two main theorems 

The first main result gives sufficient conditions on 9 for unconditional stability of each of the Do, 
CS, MCS, HV schemes in the application to for k = 2 or k = 3 under condition (|I.2[) with 
arbitrary given 7. 

Theorem 2.3 Consider equation for k = 2 or fc = 3 mt/i symmetric positive semidefinite 

matrix D and periodic boundary condition. Let 7 G [0, 1] and assume H1.2\) holds. Let il.3\) . jl.4]) 
be obtained by central second-order FD discretization and splitting as described in Section^ Then 
for the following parameter values 6 the Do, CS, MCS, HV schemes are unconditionally stable 
when applied to (1.3\) . {!■$ '■ 

• Do scheme: 

6>^{ifk = 2) and 6 > max j I, 2(27 Q + 1} | {if k = 3); 
• CS scheme: 

Q>- {if k = 2 or k = 3); 

• MCS scheme: 

9>m&x\], lill (i/ k = 2) and 9 > max (-, ^ + (i/ fc = 3); 
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• -f/V scheme: 



6>>max{-, 7 + V l(»/fc = 2) and 6>>max<[-, 27+ * 1 (ifk = 3). 

l4'4 + 2V2j W 7 " l4'4 + 2V3j W ' 

Proof The proofs for the four schemes are similar. In view of this, we shall consider here the HV 
scheme and leave the proofs for the Do, CS and MCS schemes to the reader^ 

Using the properties (|2.2h )- (|2.2b ) of the scaled eigenvalues, it is readily seen that for M = T 
the stability requirement (I1.10[) is equivalent to 

2p 2 + (2p - l){z + z) + \{z + z) 2 > , (2.5a) 
2p-l + i(z + z) >0. (2.5b) 

Condition (|2.5a ) is always fulfilled since the discriminant (2p — l) 2 — 4p 2 = — 4p + 1 < 0. Subse- 
quently, using (|2.3p and (|2.4p it is easily seen that if there exists k > such that 

k I k \ 

2l[{l + y])-l-K £^ + 7£w»i >0 (2.6) 
3=1 \3 1 »/j / 



2 See, however, [H] for the MCS scheme if k = 2. 
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for all yj > (1 < j < k), then condition (|2.5b ) is fulfilled whenever 9 > 1/(2k). 
For k = 2 the inequality (|2.6I) reads 

1 + 2{yl + yl + yivl) ~ K{y\ + y\ + 2 7 y l2/2 ) > 0, 
and this can be rewritten as 

(2 - K )( yi - y 2 ) 2 + 2 (y m + 1 - \k - i K7 ) 2 + 1 - 2 (1 - \k - i« 7 ) 2 > . 
Hence, the inequality is satisfied if 



n<2 and V2 | 1 - \k - ±kj\ < 1 



which is equivalent to 



2-V2 . J 2 + v 7 ^ 

< k < mm < 2, 



7+1 ~ ~ 1 7+1 J 
Selecting the rightmost value for k, yields that (|2.5b ) holds for k = 2 whenever 



> max 



I 7 + 1 ] 
4'4 + 2V2j 



For k — 3 the inequality (|2.6|) reads 

3 

i + (2 - K ) Vj + 2 vh) + lylvlvl - 2^7 ViVj > 0, 

j = l i<j i<j 

which, by the identity 

3 

is equivalent to 

1 + |(2 - «) ]T( yi - y,f + 2 ]T y 2 y 2 + 27/^^2 + ( 2 _ K _ 2k7 ) £ > 0. 

Let it = J/1J/2, u = 2/1J/3, w = J/22/3- Then u, u, u; > and the above inequality can be written as 

i(2-«)^( yi - yj ) 2 + 2P(u,«,w) > 0, 

i<j 

where 

P(u, v,w) = a + ?i 2 + v 2 + u) 2 + mw — 5(u + V + w) 

with 

a = i and 5 = k K + k 7 — 1 . 

Hence, if k < 2 and P(w, u, w) > 0, then (gS]) is satisfied for k = 3. If (5 < 0, i.e. rc < 2/ (27 + 1), 
then obviously P(u,v,w) > a > 0. Next consider <5 > 0. Lemma [2.11 yields that P(u,v,w) > 
whenever 

(<5 + l)(3-2v / *TT) > \ and <5 < 1. (2.7) 



Set x = V<5 + 1. Then x > 1 and 

( ( 5+l)(3-2v / 5TT) > 5 ^ 4a; 3 -6a; 2 + l < (2x-l)(2a; 2 -2a;-l) < <=^> x < 5 + 5V3. 
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It follows that (|2.7p is equivalent to S < |v3) which means k < (2 + \/3)/(27 + 1). Hence, if 



K < min 



( 




then (|2.6p holds for fc = 3. Selecting the upper bound for k, yields that (|2.5b ) is fulfilled for k = 3 
whenever 



Upon setting 7 = 1 in Theorem 12.31 the resulting sufficient conditions on 9 for the CS, MCS, 
HV schemes agree with those in [IT] Thms. 2.2, 2.5, 2.8]. The above theorem thus forms a proper 
generalization of results in [IT] . For the Do scheme, the obtained sufficient condition generalizes 
and improves the corresponding result for this scheme from [3]. In particular, if 7 = 1 and k = 3, 
then [3] yields 9 > 3^ - § w 0.696, whereas Theorem O yields 9 > §. 

In view of Theorem l2.3[ a smaller parameter value 9 can be chosen while retaining unconditional 
stability if it is known that (jl.2l) holds with certain given 7 < 1 . This is useful, in particular since 
a smaller value 9 often yields a smaller, i.e., more favorable, error constant. 

The MCS scheme with the lower bound for 9 given by Theorem 12.31 has been successfully used 
recently in the actual application to the three-dimensional Heston-Hull- White PDE from financial 
mathematics, see Haentjens & In 't Hout [5]. 

The following theorem provides, for each ADI scheme, a necessary condition on 9 for uncondi- 
tional stability. 

Theorem 2.4 Let k > 2 and 7 € [0, 1] be given and consider any ADI scheme from Section [IJ 
Suppose that the scheme is unconditionally stable whenever it is applied to a system A1.3\) , ljl-4V 
that is obtained by central second-order FD discretization and splitting as described in Section[]]of 
any equation 1 1.1)) with positive semidefinite kxk matrix D satisfying ll.S)) and periodic boundary 
condition. Then 9 must satisfy the following bound: 

• Do scheme: 




9 > max 




with dk 




• CS scheme: 



9 > max 




• MCS scheme: 



9 > max 




with bk 




I 



• HV scheme: 
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Proof As in the foregoing proof, we consider here the HV scheme and leave the (analogous) proofs 
for the other three ADI schemes to the reader. 

Consider the k x k matrix D = (<%) with da = 1 and = 7 whenever 1 < i ^ j < k. 
Clearly, D is symmetric positive semidefinite and (|1.2p holds. Let Axi = Ax > 0, so that 
m =r:= At/(Ax) 2 . 

First, choose the angles 4>j in (|l-17j> equal to zero for j > 2 . Then the eigenvalues Zj are given 
by zo = 0, z\ = — 2r(l — cos0i) and z-i — Z3 = . . . = Zk = 0. In view of (|2.5b ). we have 

2p-\ + \(zo + z) = \ + (\- 26)z x > 

whenever r > 0, <f>x € [0, 2w). This immediately implies 9 > j. 

Next, choose all angles <j>j in (jl . 1 7[) the same, i.e., (j>j = (j). Then the eigenvalues Zj are given 

by 

z = -rk(k - l)7[sin 2 4> - - cos0) 2 ] , 
Zj = -2r(l-cos0), i = l,2,...,fc, 



where 



fe(fc- 1) ■ 

Note that |/3| < 1. By P3b). it must hold that 

10z o + 6»z fc (k - l) 7 6»r[ sin 2 <f>-~P(l - cos^) 2 ] + 26»r(l - cose 
9 > 



2 2p-l 2 2[l + 26»r(l-cos0)] fe - 1 

whenever r > 0, G [0, 27r). Setting a = 29r{\ — cos0), this yields 

fe (fc - l)7a[l + cos^-^(l - cos0)]/2 + a 
~ 2 2(1 + a) k - 1 

for all a > 0, <p € (0, 27r). Taking the supremum over cj), gives 

9 > - h(a) [(k - 1)7 + 1] for all a > 0, 



2 



where 



h(a) 



ka 



2(l + a) fe - 1 



By elementary arguments (cf. [TT]) it follows that the function h has an absolute maximum which 
is equal to given in the theorem. 



Upon setting 7 = 1 in Theorem 12.41 the necessary conditions on 9 for the CS, MCS, HV 
schemes reduce to those given in [11] Thms. 2.3, 2.6, 2.9]. Further, for the Do scheme and 7 = 1 
there is agreement with the necessary condition from [3]. 

It is readily verified that, for each ADI scheme, the sufficient conditions of Theorem 12.31 and 
the necessary conditions of Theorem 12.41 are identical whenever k = 2 or k — 3 and < 7 < 1 . 
Hence, in two and three spatial dimensions, these conditions are sharp. 

The interesting question arises whether the necessary conditions of Theorem 12.41 are also suf- 
ficient in spatial dimensions k > 4. In |llj it was proved that this is true for the HV scheme 
and 7 = 1. It can be seen, however, that the proof from |11) does not admit a straightforward 
extension to values 7 < 1. Further, in the case of the Do, CS, MCS schemes a proof is not clear 
at present. Accordingly, we leave this question for future research. 
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3 Numerical illustration 



In this section we illustrate the main results of this paper, Theorems 12.31 and 12.41 We present 
experiments where all four ADI schemes (|1.5p - (|1.8p are applied in the numerical solution of mul- 
tidimensional diffusion equations possessing mixed derivative terms. The PDE is semidis- 
cretized by central second-order finite differences as described in Subsection II .31 with = 0, and 



the scmidiscrctc matrix A is splittcd as described in Subsection 11.21 The boundary condition is 
taken to be periodic (in this case g(t) = 0). 

The first experiment deals with (|1.1[) in two spatial dimensions. We choose initial function 

U(X 1 ,X 2 ,0) = g-^in^iHcos 2 ^)] (0 < XuX2 < X ) 



and diffusion matrix D given by 



D = 0.025 



1 2 7 
2 7 4 



This matrix D is positive semidefinite and the condition (|1.2[) holds whenever 7 6 [0,1]. Consider 
7 = 0.9. Then the lower bounds on 9 provided by Theorem 12.31 for the MCS and HV schemes 
are, rounded to three decimal places, equal to 0.317 and 0.278, respectively. For the Do and CS 
schemes, the lower bound always equals h in two spatial dimensions, independently of 7. 

Figure Q] shows for Ax\ = Ax2 = 1/80 the semidiscrete solution values U (0) and U(5) displayed 
on the grid in f2, so that they represent the exact solution u at t = and t = 5. 



t = 



t = 5 




0.095 




Figure 1: Exact solution u of 2D problem at t = and t = 5. Note the different scales on the 
vertical axes. 



To the semidiscrete systems with Ax± = Ai2 = 1/m for m = 40, 80 we have applied the Do, 
CS, MCS, HV schemes using their exact lower bound values 9 (above) as well as the values 0.45, 
0.45, 0.29, 0.25, respectively, which are approximately 90% of these. For a sequence of step sizes 
At = 1/N with 10 -3 < At < 10° we computed the global temporal errors at time t = 5, defined by 

e(At;m) = m- 1 \U(5)-U 5N \ 2 , 

where | ■ 1 2 is the Euclidean vector norm and 1/m is a normalization factor. Figure [2] displays the 
obtained result. 

From the right column of Figure [5] it is clear that, for each ADI scheme, using the smaller 
value 9 leads to large temporal errors for natural step sizes. We verified that these large errors 
correspond to the spectral radii of the pertinent iteration matrices being greater than 1. This 
indicates that they are caused by a lack of stability. Additional experiments with larger values t 
(e.g. t = 10) also shows that the large temporal errors are amplified, as one may expect. 
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m = 40 m = 40 




Figure 2: Global temporal errors e(Ai; m) versus At for Do, CS, MCS, HV schemes when applied 
to a semidiscretized 2D problem (II. ip with 7 = 0.9. Top row: m — 40. Bottom row: m = 80. 
Left column: lower bound values 9 given by Theorem [273] Right column: values 9 that are about 
90% of these. 



The left column of Figure [2] displays the results in the case where the lower bound values 9, 
given by Theorem l2.31 are used. Then all temporal errors are bounded from above by a moderate 
constant and decrease monotonically when At decreases. This suggests an unconditionally stable 
behavior of the schemes. A further examination in this case shows that the global temporal errors 
for the Do scheme can be bounded from above by CAt and for the CS, MCS, HV schemes by 
C(At) 2 (whenever At > 0) with constants C depending on the scheme. This clearly agrees with 
the respective orders of consistency of the schemes. Moreover, for each scheme the constant C is 
only weakly dependent on the number of spatial grid points, determined by m, indicating that the 
error bounds are valid in a stiff, hence favorable, sense. We remark that similar conclusions were 
found for larger values t. Further we numerically verified that upon increasing 9 from its lower 
bound value, the error constant C increases, as mentioned in the discussion following Theorem 12. 31 

The second experiment deals with Ijl.ljl in three spatial dimensions. We take initial function 

u(xi,a;a,a;s,0) = e-^^O+cos^^WV^)] ( < x u x 2 ,x 3 < 1) 
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m = 40 m = 40 




Figure 3: Global temporal errors e(Af; m) versus At for Do, CS, MCS, HV schemes when applied 
to a semidiscretized 3D problem (jl.ll) with 7 = 0.75. Top row: m = 40. Bottom row: m = 80. 
Left column: lower bound values 9 given by Theorem [273] Right column: values 9 that are about 
90% of these. 



and diffusion matrix D given by 

/ 1 2 7 7 \ 
D = 0.025 2 7 4 27 , 

V 7 2 7 1 j 

which is positive semidefinite and such that condition (|1.2[) holds whenever 7 <E [0, 1]. Here we 
take 7 = 0.75. Then the lower bounds on 9 given by Theorem 12.31 for the Do, MCS, HV schemes 
are, rounded to three decimal places, equal to 0.556, 0.385, 0.335, respectively. For the CS scheme 
the lower bound is always \ in three spatial dimensions, independently of 7. To the semidiscrete 
systems with Axi = Ax-i = Ax 3 = 1/m for m = 40, 80 we have applied the Do, CS, MCS, HV 
schemes using their exact lower bound values 9 as well as the values 0.5, 0.45, 0.35, 0.3, respectively, 
which are approximately 90% of these. Figure [3] displays the normalized global temporal errors 

e(Ai;m) = m~ 3/2 |C/(5) - U 5N \ 2 

for a sequence of step sizes At — l/N with 10 -3 < At < 10°. Exactly the same observations can 
be made as in the experiment (above) for the two-dimensional case. The large temporal errors for 
each ADI scheme when applied with the smaller value 9, as seen in the right column of Figure [3] 
correspond to instability of the scheme. When applied with their lower bound values 9, given by 



13 



Theorem 12. 31 the error behavior for all ADI schemes is in agreement with unconditional stability 
of the schemes. Moreover, in this case a stiff order of convergence is observed that is equal to one 
for the Do scheme and equal to two for the CS, MCS and HV schemes. 

4 Conclusions 

In this paper we analyzed stability in the von Neumann sense of four well-known ADI schemes - the 
Do, CS, MCS and HV schemes - in the application to multidimensional diffusion equations with 
mixed derivative terms. Such equations are important, notably, to the field of financial mathemat- 
ics. Necessary and sufficient conditions have been derived on the parameter for unconditional 
stability of each ADI scheme by taking into account the actual size of the mixed derivative terms, 
measured by the quantity 7 G [0, 1]. Our two main theorems generalize stability results obtained 
by Craig & Sneyd [3] and In 't Hout & Welfert [TT]. Ample numerical experiments have been 
presented, illustrating the main theorems and also providing insight into the convergence behavior 
of all schemes. 

Among issues for future research, it is of much interest to derive sufficient conditions on 9 for 
unconditional stability of each ADI scheme in arbitrary spatial dimensions k > 4 and arbitrary 
7 £ [0, 1]. Also, it is of much interest to extend the results obtained in this paper to equations 
with advection terms. This leads to general complex, instead of real, eigenvalues, which forms a 
nontrivial feature for the analysis, cf. e.g. [51 [TP]. 
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